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Abstract.  In  the  plumes  of  Hall  thrusters  and  ion  thrusters,  high  energy  ions  experience  elastic  collisions  with  slow  neutral  atoms. 
These  collisions  involve  a  process  of  momentum  exchange,  altering  the  initial  velocity  vectors  of  the  collision  pair.  In  addition  to 
the  momentum  exchange  process,  ions  and  atoms  can  exchange  an  electron,  resulting  in  slow  charge-exchange  ions  and  fast  atoms. 
In  these  simulations,  it  is  particularly  important  to  accurately  perform  computations  of  ion-atom  elastic  collisions  in  determining 
the  plume  current  profile  and  assessing  the  integration  of  spacecraft  components.  The  existing  models  are  currently  capable  of 
accurate  calculation  but  are  not  fast  enough  such  that  the  calculation  can  be  a  bottleneck  of  plume  simulations. 

This  study  investigates  methods  to  accelerate  an  ion-atom  elastic  collision  calculation  that  includes  both  momentum-  and 
charge-exchange  processes.  The  scattering  angles  are  pre-computed  through  a  classical  approach  with  ab  initio  spin-orbit  free 
potential  and  are  stored  in  a  two-dimensional  array  as  functions  of  impact  parameter  and  energy.  When  performing  a  collision 
calculation  for  an  ion-atom  pair,  the  scattering  angle  is  computed  by  a  table  lookup  and  multiple  linear  interpolations,  given  the 
relative  energy  and  randomly  determined  impact  parameter.  In  order  to  further  accelerate  the  calculations,  the  number  of  collision 
calculations  is  reduced  by  properly  defining  an  effective  elastic  collision  cross-section  that  is  used  to  decide  if  the  degree  of 
momentum  exchange  is  insignificant.  In  the  MCC  method,  the  target  atom  needs  to  be  sampled;  however,  it  is  confirmed  that  initial 
target  atom  velocity  does  not  play  significant  role  in  typical  electric  propulsion  plume  simulations  such  that  the  sampling  process 
is  unnecessary.  With  these  implementations,  the  computational  run-time  to  perform  a  collision  calculation  is  reduced  significantly 
compared  to  previous  methods,  while  retaining  accuracy  of  the  high  fidelity  models. 


INTRODUCTION 

In  the  plumes  of  Hall  effect  thrusters  and  ion  thrusters,  ions  electrostatically  accelerated  to  high  energy  can  experience 
elastic  collisions  with  the  slow  neutral  atoms  exiting  the  thrusters.  During  the  collisional  process,  the  high  energy 
ions  exchange  momentum  with  their  collision  partners,  and  some  fraction  of  the  population  is  deflected  at  angles 
greater  than  the  plume  divergence.  Meanwhile,  the  elastic  collision  may  involve  an  exchange  of  one  or  more  electrons, 
leading  to  a  slow  CEX  ion  and  a  high  energy  atom  after  each  charge-exchange  (CEX)  process.  Therefore,  the  CEX 
collision  is  considered  as  a  subset  of  elastic  collision,  and  we  call  the  elastic  collision  without  an  exchange  of  charge 
a  MEX  collision.  Proper  modeling  of  these  collisions  is  particularly  important  in  accurately  determining  the  plume 
current  profile  and  assessing  the  integration  of  spacecraft  components  as  well  as  the  thruster  life  and  long  duration 
performance;  the  slow  CEX  ions  are  prone  to  electric  fields  that  direct  the  ions  toward  the  spacecraft  and  thruster 
components,  resulting  in  contamination  through  sputtering  and  deposition  after  gaining  significant  energies. 

The  elastic  collision  between  high  energy  ions  and  slow  atoms  has  been  approximated  by  a  few  different  methods. 
In  a  particle-in-cell  (PIC)  model  developed  by  Oh  [1],  simple  models  were  employed;  the  elastic  collision  was  treated 
as  isotropic  scattering  whereas  the  CEX  collision  was  handled  without  any  exchange  of  momentum.  More  recently, 
the  elastic  collision  has  been  calculated  by  a  detailed  model  that  solves  the  classical  scattering  equation  with  a  precise 
interatomic  potential  [2,  3,  4].  This  method  provides  much  more  accurate  scattering  characteristics  necessary  for 
plume  simulations.  However,  the  numerical  integration  of  classical  scattering  equation  is  computationally  expensive 
such  that  the  collision  calculation  can  be  a  bottleneck  to  a  massive  particle  simulation.  An  alternative  method  is  to 
use  a  curve-fit  representation  of  a  center-of-mass  (CM)  differential  cross-section  to  determine  post-collision  particle 
velocities  [5,  6].  This  method  is  much  faster  than  the  direct  numerical  integration.  However,  it  is  limited  to  the  energy 
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that  the  differential  cross-section  is  provided  at.  At  the  exit  plane  of  the  Hall  thruster,  ions  do  not  necessarily  have  the 
energy  equivalent  to  the  potential  difference  between  a  cathode  and  an  anode,  as  they  are  created  at  different  locations 
along  the  channel.  Furthermore,  their  energies  after  the  first  collision  are  always  reduced  significantly  for  the  case  of 
large  angle  scattering.  Therefore,  the  collision  calculation  method  does  not  provide  accurate  collision  characteristics 
for  every  particle  within  the  plume. 

The  objective  of  this  research  is  to  develop  a  fast  and  accurate  method  to  compute  the  scattering  angle  for  an 
elastic  collision  between  high  energy  ions  and  slow  atoms.  Instead  of  using  the  experimentally  determined  differential 
cross-section,  the  CM  differential  cross-section  can  be  numerically  computed  from  the  classical  theory  at  any  incident 
ion  energy.  By  determining  an  applicable  range  of  impact  parameters  for  a  given  energy  and  formulating  a  two- 
dimensional  array  of  scattering  angles  as  functions  of  impact  parameter  and  energy,  the  scattering  angle  can  be  quickly 
obtained  by  a  table  lookup.  This  paper  provides  the  detail  of  the  table  formulation  and  compares  its  performance  with 
other  methods.  Although  the  paper  deals  primarily  with  Xe+  +  Xe  collision,  the  same  concept  can  be  applied  to 
collisions  between  other  species. 


NUMERICAL  METHOD  FOR  COLLISION  CALCULATION 


Classical  Scattering  Model 

Two  classes  of  collision  models,  Monte  Carlo  Collision  (MCC)  and  Direct  Simulation  Monte  Carlo  (DSMC)  models, 
are  commonly  used  for  particle  simulations  of  plasma  with  collisions.  Both  models  involve  determination  of  scattering 
angles  and  alteration  of  particle  velocity  vector  when  applying  elastic  collisions  to  particles.  Phenomenological  models 
such  as  variable  hard-sphere  (VHS)  model  [7]  have  been  successfully  used  in  simulations  of  rarefied  gas.  These  models 
are  designed  to  reproduce  the  property  of  real  gases  by  choosing  adjustable  parameters  to  match  the  experimental  data. 
Furthermore,  simple  scattering  characteristics  allow  fast  computation  of  post-collision  velocity  vectors.  In  plumes  of 
spacecraft  propulsion  devices,  the  background  gas  density  rapidly  reduces  away  from  the  thrusters  so  that  at  most  a 
few  elastic  collision  events  per  incident  particle  would  take  place.  Therefore,  the  unrealistic  scattering  characteristics 
employed  in  the  phenomenological  models  would  not  permit  in  accurate  determination  of  the  plume  current  profile. 
The  scattering  characteristics  of  ions  and  neutral  atoms  can  be  obtained  more  precisely  through  the  classical  approach, 
integrating  the  interatomic  force  between  the  collision  pair  along  the  incident  particle’s  trajectory  as  shown  in  Figure  1. 
The  deflection  angle  in  the  center  of  mass  (CM)  frame,  x*  is  expressed  as 


X(b,  Er)  =  n-  2b 


f 


dr 


Rm(b,Er )  r2[\-b2/r2-V(r)/Er]2 


(1) 


where  b  is  the  impact  parameter,  r  is  the  internuclear  distance,  V(r)  is  the  interatomic  potential,  and  Er  is  the  relative 
energy.  The  classical  turning  point  or  the  minimum  distance  of  approach,  Rm ,  is  calculated  by  finding  the  largest  root 
of  the  equation. 

1  -  b2/Rl  -  V(Rm)/Er  =  0  (2) 

Equation  (1)  and  (2)  may  be  evaluated  analytically  if  the  potential  function  V(r)  is  simple.  However,  this  is  often 
not  the  case  for  a  function  that  closely  approximates  the  real  interaction  potential.  Computing  Equation  (1)  with 
direct  numerical  integration  can  be  difficult  since  the  denominator  approaches  zero  in  the  limit  of  r  =  Rm.  Instead, 
Equation  (1)  is  solved  by  applying  the  Gauss-Mehler  formula  [8]. 


FIGURE  1.  Classical  scattering  trajectory. 
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FIGURE  2.  Spin-orbit  free  interaction  potential  energy  curve  for  Xe+(2P)+Xe  calculated  by  Paidarova  and 
Gadea  [10]. 


The  interaction  potentials  for  Xe+  +  Xe  collisions  have  been  calculated  with  ab  initio  quantum  chemistry  [9, 
10,  11].  In  the  computational  model,  four  spin-orbit  free  potentials,  11  and  £  potentials  with  the  gerade  (g)  and  the 
ungerade  (u)  states,  calculated  by  Paidarova  and  Gadea  are  used  [10].  These  potentials  are  shown  in  Figure  2.  Chiu 
et  al.  [12]  verified  that  the  scattering  results  are  not  affected  significantly  by  using  the  full  spin-orbit  potentials.  The 
averaged  (u,g)  pairs  of  the  spin-orbit  free  potentials  are  fitted  using  Morse  potential  form  and  the  fitting  parameters 
are  given  in  Ref.  12. 

Vavg(r)  =  De(e2b^~r)  -  2eb^-r))  (3) 

The  statistical  weights  for  £  and  II  potentials,  d %  and  dn,  are  1/3  and  2/3,  respectively. 

Figure  3  shows  the  center  of  mass  deflection  angle  as  a  function  of  impact  parameter  for  different  center  of  mass 
energies.  Note  that  Er-  150  eV  roughly  corresponds  to  the  laboratory  (LAB)  frame  incident  ion  energy  ( E  =  300  eV) 
typically  seen  in  a  Hall-effect  thruster  plume.  The  deflection  function  for  Er  =  0.5  eV  has  a  small  minimum  at  b  «  4.9 
A,  which  gives  a  rainbow  singularity  in  the  differential  cross-section  at  the  rainbow  angle  corresponding  to  the  mini¬ 
mum  deflection  angle.  Also,  the  deflection  function  becomes  discontinuous  at  the  critical  angle  at  even  lower  energies 
(i.e.  orbiting  singularity).  These  singularities  are  not  realistic,  and  quantum  scattering  has  to  be  considered  [13,  14]. 
However,  at  high  energies,  the  deflection  is  mostly  affected  by  the  repulsive  part  of  the  potential  at  short-range  so  that 
the  deflection  function  barely  has  a  minimum.  In  this  regime,  the  classical  scattering  model  is  sufficiently  accurate 
to  approximate  ion-atom  collision  behavior  [14].  For  all  the  energies,  the  CM  scattering  angle  approaches  zero  as 
b  — >  oo,  however,  it  never  actually  reaches  zero.  This  implies  that  the  total  collision  cross-section  is  unbounded  unless 
the  quantum  scattering  theory  is  considered.  In  applying  the  classical  model,  it  is  necessary  to  provide  a  cut-off  impact 
parameter  or  deflection  angle  [7].  The  cut-off  impact  parameter  is  chosen  from  the  deflection  functions  at  different 
energies  as  discussed  hereinafter. 


Total  Cross-Section 

Charge  exchange  collision  cross-sections  have  been  measured  experimentally  and  fitted  to  a  simple  logarithmic  for¬ 
mula  by  Miller  et  al.  [15]. 

o-cex  =  87.3  -13.6  log  (E)  (4) 

Equation  (4)  generally  agrees  well  with  previous  experimental  measurements  and  has  been  verified  theoretically  by  a 
semi-classical  calculation.  The  semi-classical  calculation  involves  an  integration  of  charge-exchange  probability,  Fcex 
(see  Figure  4),  evaluated  from  the  elastic  phase  shifts  for  the  two  states.  At  small  impact  parameters,  Fcex  oscillates 
rapidly  between  0  and  1,  and  the  average  value  of  0.5  is  typically  used  in  models.  At  larger  impact  parameters, 
its  maximum  value  in  the  oscillation  starts  to  decrease  due  to  the  stronger  interaction  of  2  potentials.  However, 
the  deflection  is  close  to  zero  in  this  regime  such  that  precise  approximation  of  charge  exchange  probability  is  not 
necessary.  For  this  reason,  previous  computational  models  have  approximated  elastic  collision  cross-section  to  be 
twice  the  CEX  cross-section,  i.e.  <xtotai  =  2<xcex,  and  used  Fcex  =  0.5  [2,  3]. 


FIGURE  3.  Deflection  functions  computed  at  different 
CM  energies. 
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FIGURE  4.  Charge-exchange  probability  for  a  relative 
energy  of  150  eV. 


Larger  cross-section  leads  to  more  collision  calculations  due  to  the  increased  collision  probability.  Calculation 
of  post-collision  velocity  involves  computation  of  scattering  angle  and  multiple  coordinate  transformations  between 
LAB  and  CM  frames  and  between  pre-  to  post-collision  velocity  reference  frames.  Performing  these  calculations  for 
a  number  of  collision  pairs  can  be  computationally  expensive.  Therefore,  the  computational  cost  can  be  reduced  by 
simply  skipping  the  collision  calculation  for  collision  pairs  that  do  not  experience  significant  momentum  transfer. 
Let  us  define  an  effective  elastic  collision  cross-section,  <xe i,  which  acts  as  a  limit  that  deflection  can  occur.  Note 
that  the  CEX  process  may  still  occur  when  the  impact  parameter  is  within  <xei.  The  cross-section  can  be  obtained  by 
choosing  the  cut-off  impact  parameter  from  the  deflection  functions  at  different  energies.  Finding  the  impact  parameter 
corresponding  to  the  cut-off  deflection  angle  of  1°  and  calculating  the  corresponding  cross-section  by  cr  =  nb 2,  <xei 
can  be  fitted  to  a  logarithmic  formula. 

cr ei  =62.3  -13.4  log  (E)  (5) 

The  cut-off  deflection  angle  can  be  chosen  to  be  other  than  1°,  but  a  larger  angle  would  result  in  a  loss  of  accuracy 
while  reducing  the  computational  cost.  The  classical  theory  is  no  longer  accurate  within  the  glory  singularity  regime 
such  that  a  smaller  cut-off  angle  would  not  necessarily  increase  the  accuracy.  At  the  LAB  energy  of  300  eV,  <xtotai  = 
107.2  A  and  <xei  =  29.1  A;  therefore,  roughly  73%  of  collision  calculation  can  be  cut. 


Differential  Cross-Section 

Although  the  differential  cross-section  is  not  used  anywhere  in  the  collision  calculation,  it  is  covered  herein  to  facilitate 
later  discussions.  The  differential  cross-section  is  useful  in  verifying  the  implementation  as  it  can  be  obtained  by 
(a)  numerically  evaluating  from  theory  and  (b)  applying  collisions  and  sampling  particles  in  the  code.  Unlike  the 
deflection  function  shown  in  Figure  3,  the  differential  cross-section  includes  both  the  MEX  and  CEX  collision  effect. 
Once  the  deflection  angles  at  different  impact  parameters  are  computed  by  Equation  (1),  the  differential  cross-section 
can  be  determined  by  the  following  equation. 


Icm  (T>  Er ) 


dcr(x ) 

b 

d£l 

CM 

sin  x(dx/db) 

(6) 


Here,  dx/db  can  be  approximated  by  using  first  order  difference  equations.  Taking  into  account  backscattered  CEX 
ions,  the  differential  cross-section  can  be  rewritten  as  [12] 


^cm(l)  =  (1  -  Pcex(x))Icm (x)  +  PcexIcm(tt  ~x) 


(7) 


The  first  and  second  terms  on  the  right-hand  side  of  Equation  (7)  denote  the  MEX  and  CEX  collision  contributions 
to  the  total  differential  cross-section,  respectively.  If  the  collision  partners  have  equal  mass  and  the  target  particle  is 


stationary,  the  CM  differential  cross  sections  can  be  converted  into  the  LAB  reference  frame  by  [2] 


dcr 


LAB 


_  ^CM  (L)4  cos(y/2) 


(8) 


For  the  CEX  collision,  the  LAB  differential  cross  section  is  calculated  by  replacing^  with  n  ~x  in  Equation  (8).  The 
high  probability  in  small  angle  scattering  for  the  energy  of  interest  leads  to  large  differential  cross-sections  at  extreme 
angles  of  0°  and  90°  as  a  result  MEX  and  CEX  collisions,  respectively. 


Formulation  of  Matrix  for  Table  Lookup 

In  determining  the  scattering  angle  by  computing  Equation  (1),  the  Gauss-Mehler  formula  works  significantly  better 
than  the  direct  numerical  integration  with  respect  to  speed  and  accuracy.  However,  performing  the  calculation  for 
a  number  of  particles  on  the  fly  can  be  computationally  expensive.  A  much  more  attractive  way  is  to  pre-compute 
the  scattering  angles  for  ranges  of  energies  and  impact  parameters  and  perform  a  linear  interpolation  of  the  stored 
scattering  angles.  In  the  2-dimensional  array  of  scattering  angles,  the  energy  values  are  distributed  uniformly  within 
the  minimum  and  maximum  CM  energies,  defined  as  Emm  and  Ema x,  respectively. 

Ej  =  Emin  +  j  A£,  0  <  j  <  Ne  -  1  (9) 

where  A E  =  (Emax  -  Emin)/(NE  -  1)  and  Ne  is  the  number  of  energy  values.  Each  discrete  energy  value,  Ej,  is 
associated  with  a  cut-off  impact  parameter,  bmax ,  defined  by  the  effective  elastic  cross-section  given  in  Equation  (5). 
The  minimum  impact  parameter  is  simply  zero  so  that  the  impact  parameter  values  to  evaluate  the  scattering  angle  are 
uniformly  distributed  between  0  and  bmax. 

bi  =  i  bmax/(Nb  -  1),  0  <  i  <  Nb  -  1  (10) 

The  scattering  angle,  xuj*  is  evaluated  for  every  combination  of  bt  and  Ej  and  is  stored  in  the  2-dimensional  array.  In 
order  to  ensure  a  smooth  variation  of  scattering  angle  near  the  cut-off  impact  parameter,^  is  forced  to  be  zero  at  bmax. 

For  a  given  pair  of  particles,  the  CM  energy  is  not  necessarily  at  the  values  used  to  formulate  a  table.  Therefore, 
the  scattering  angle  is  determined  by  interpolating  to  the  values  for  the  collision  pair.  First,  the  fractional  energy  index, 
j,  for  the  table  is  obtained  such  that  the  CM  energy,  Ep ,  is  bounded  by  Ej  and  Ej+\. 


jf  =  (Ep-Emin)/AE  (11) 

The  integer  value  of  jf  is  used  as  the  energy  index,  i.e.  j  =  int (jf),  and  the  weight  for  linear  interpolation  is  simply 
given  as  Wj  =  jf  -  j .  Since  the  maximum  impact  parameter  differs  for  every  energy,  the  fractional  impact  parameter 
indices  are  not  the  same  for  Ej  and  Ej+i- 

ifj  =  bp/Abj  (12) 
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FIGURE  5.  Illustration  of  the  interpolation  process  used  to  determine  scattering  angle.  Indices  of  m  and  n  are  used  in 
place  of  ij  and  ij+ 1,  respectively. 


(a)  (b) 

FIGURE  6.  LAB  differential  cross-section  computed  from  theory  and  by  sampling  particles  from  the  code.  The 
impact  parameter  values  in  the  table  is  distributed  (a)  uniformly  in  b  and  (b)  uniformly  in  b2. 


Indices  and  weights  for  ij  and  ij+\  are  obtained  similarly  as  the  case  of  energy.  Finally,  the  scattering  angle  is  obtained 
by  performing  three  linear  interpolations. 

X(bp,  Ej)  =  (1  -  WiJxibij,  Ej)  +  WijX(bij+uEj) 

X(bp ,  Ej+ 1)  =  (1  -  wij+l)x(bij+1,Ej+  0  +  wij+lx(bij+l+uEj+1)  (13) 

X(bp,  Ep)  =  (1  -  Wj)x(bp,  Ej)  +  wjxibp,  Ej+1) 

The  interpolation  process  is  also  illustrated  in  Figure  5. 

Tabulation  of  scattering  angle  and  use  of  a  look-up  table  have  also  been  implemented  by  Sharipov  and  Strapas- 
son  [16].  However,  the  distribution  of  impact  parameter  is  different  in  that  they  used  uniform  distribution  of  b2  instead 
of  b.  Their  approach  ensures  equal  probability  that  the  impact  parameter  lie  within  every  bin  (e.g.  between  bj  and 
bj+ 1).  Unlike  their  application,  the  most  important  collision  within  the  spacecraft  plume  is  at  high  energy  such  that  the 
deflection  function  is  roughly  monotonically  decreasing  with  larger  impact  parameters  (see  Figure  3).  In  order  to  cap¬ 
ture  the  plume  current  profile  correctly,  it  is  necessary  to  resolve  the  small  impact  parameter  region  where  large  angle 
scattering  occurs.  Due  to  the  steeper  slope  in  deflection  function  as  b  approaches  b  -  0  A,  a  uniform  distribution  in  b2 
results  in  less  accurate  scattering  angles  near  the  region.  Therefore,  a  uniform  distribution  in  b  is  chosen  to  improve 
the  accuracy  as  b  approaches  b  =  0  A  for  a  given  number  of  impact  parameter  in  the  table.  Figure  6  compares  the 
differential  cross-section  computed  from  the  code  by  sampling  particles  after  a  single  elastic  collision.  In  constructing 
a  table  to  determine  scattering  angles,  50  impact  parameter  values  are  used  (Nb  =  50).  It  is  readily  seen  that  the  uni¬ 
form  distribution  of  b  agrees  much  better  with  the  theory  than  a  uniform  distribution  of  b2.  The  sampled  differential 
cross-section  approaches  the  theoretical  one  as  Nb  is  increased. 

Implementation  in  Monte  Carlo  Collision  Method 

Figure  7  shows  the  flowchart  of  the  elastic  collision  calculation  algorithm,  utilizing  three  different  cross-sections  (i.e. 
cr total  ?  ctq i,  and  <xcex)  given  in  Equation  (4)  and  (5).  The  total  cross-section,  <xt0 tai,  can  be  arbitrarily  chosen  as  long 
as  it  is  greater  than  <xcEX+0.5<xei.  The  case  shown  in  Figure  7  assumes  that  <xtotai  =  2<xcex  such  that  Pcex  =  0.5 
for  any  value  of  impact  parameter.  For  a  case  where  <xtotai  ^  2<xcex,  Fcex  for  b  >  bmax  needs  to  be  refined  such 
that  cr cex  =  FcEX,icrei  +  Fcex,2  (crtotai  -  cre i)  is  maintained,  where  Pcex,i  =  0.5  covers  b  <  bmax  and  Fcex,2  is  the 
adjustable  parameter  that  covers  b  >  bmax. 

The  Monte  Carlo  collision  routine  involves  a  step  to  determine  whether  the  collision  between  the  collision  pair 
occurs  by  comparing  a  collision  probability,  Ptotai ,  with  a  random  number,  U \ ,  between  0  and  1 , 


P total  =  1  -  exp(-y^oCr total Vrg/Af) 


(14) 


FIGURE  7.  Simplified  flowchart  of  elastic  collision  calculation  algorithm  in  context  of  MCC  collision  method. 


where  y  is  the  collision  factor  that  is  1  for  the  same  species  pair  and  0.5  for  unlike  species  pair,  no  is  the  target  gas 
density,  vrei  is  the  relative  speed  between  the  collision  pair,  and  At  is  the  time-step.  If  Ptotai  >  U\  is  satisfied,  the 
next  step  involves  determining  if  the  scattering  angle  is  significant  enough  to  be  computed.  The  impact  parameter  for 
the  collision  event  is  b  -  bmax  ^[Th  where  U2  is  another  random  number.  If  the  impact  parameter  is  smaller  than  the 
radius  of  effective  elastic  collision  cross-section  (i.e.  b  <  'crei  /tt ),  then  the  post-collision  velocity  is  computed  using 
the  impact  parameter.  The  post  collision  velocity  calculation  is  done  in  the  CM  frame  so  that  it  involves  multiple 
coordinate  transformations  and  a  scattering  angle  calculation  by  a  table  lookup  as  discussed  in  the  previous  sections. 
After  the  elastic  collision  calculation,  the  velocity  vector  is  switched  if  the  third  random  number,  C/3,  is  smaller  than 
the  charge  exchange  probability. 

The  overall  cost  of  the  collision  calculation  is  minimized  if  the  smallest  possible  total  cross-section  is  used,  as 
it  leads  to  the  smallest  value  for  Ptotai  according  to  Equation  (14)  such  that  the  chance  of  experiencing  a  collision  is 
minimized.  In  this  study,  <xtotai  =  2<xcex  is  chosen  instead  of  <xcEX+0.5<xei  due  to  the  simplicity  in  implementation, 
but  the  refinement  can  be  done  easily  with  a  slight  modification. 


COMPARISON  OF  PERFORMANCE 

The  collision  calculation  methods  are  compared  in  context  of  their  performance  in  speed.  In  this  study,  approximately 
1, 146, 700  particles  are  injected  into  the  domain  at  a  speed  of  19, 920  m/s  (270  eV).  The  target  gas  density  is  set  to 
a  very  high  value  such  that  all  the  particles  undergo  one  elastic  collision  per  step  (i.e.  Ptotai  ~  1).  The  calculation  is 
repeated  100  times,  and  only  the  time  to  perform  the  collision  calculation  is  reported.  We  used  the  implementation  in 
COLISEUM  [5]  as  a  baseline  and  incrementally  added  improvements  including  (1)  a  lookup  table  and  (2)  an  effective 
elastic  collision  cross-section.  Lastly,  we  performed  the  calculation  assuming  (3)  a  stationary  background  target  gas, 
after  implementing  the  two  improvements.  All  of  these  implementations  were  checked  to  assess  their  correctness  by 
sampling  the  particle’s  scattering  angle  and  comparing  with  the  differential  cross-section.  The  current  model  with 
the  MCC  method  samples  the  target  particle  velocity  assuming  the  Maxwellian  distribution.  However,  the  actual 
distribution  in  a  plume  should  be  close  to  the  one  obtained  by  tracing  neutral  particles  leaving  the  thruster  exit  face 
at  a  cosine  distribution.  Nevertheless,  the  target  particle  speed  is  so  slow  that,  after  becoming  an  CEX  ion,  its  particle 


TABLE  1.  Performance  comparison  of  collision  calculation  methods.  (1),  (2),  and  (3)  are  incrementally  implemented 


into  the  code. 


Wall  Time,  s 

%  Baseline 

Baseline 

61.67 

100.0 

(1)  Lookup  Table 

45.75 

74.2 

(2)  Effective  Elastic  Collision  cr 

36.21 

58.7 

(3)  Stationary  Background  Gas  Assumption 

17.63 

28.6 

motion  is  steered  by  the  electric  field,  and  the  particle’s  initial  velocity  does  not  affect  its  trajectory  significantly.  Note 
that  this  would  not  be  an  issue  in  the  DSMC  method  as  it  does  not  require  any  sampling  from  a  fluid  property  but 
rather  requires  a  selection  of  target  particles  from  a  simulation  particle  list.  It  should  also  be  noted  that  the  assumption 
of  stationary  target  particle  only  works  for  the  case  where  the  electric  potential  within  the  plume  is  significantly  higher 
than  the  target  gas  energy,  which  is  often  the  case  for  electric  propulsion  plumes.  The  validity  of  the  stationary  gas 
assumption  in  a  plume  simulation  is  further  discussed  in  the  next  section. 

Table  1  shows  the  wall  time  for  the  collision  calculation  and  its  percentage  relative  to  the  baseline  method. 
While  the  baseline  method  is  nearly  optimized  in  sampling  the  scattering  angles  from  the  differential  cross-section, 
the  acceptance-rejection  routine  used  in  the  baseline  method  was  still  significantly  more  expensive  compared  to  a 
simple  lookup  from  a  table,  as  indicated  by  the  percent  reduction  in  wall  time  of  25.8  percent.  Furthermore,  with 
an  implementation  of  <xei,  an  additional  15.5  percent  reduction  in  wall  time  was  achieved.  At  this  point,  the  most 
expensive  operation  was  to  sample  the  target  particle  velocity  from  a  Maxwellian  distribution,  even  though  this  would 
not  lead  to  a  more  accurate  solution  in  any  way  from  the  stationary  background  gas  assumption.  By  skipping  the 
sampling  and  assuming  a  target  gas  velocity  of  zero,  another  30.1  percent  of  collision  calculation  time  was  cut.  The 
overall  reduction  of  the  collision  calculation  was  71.4  percent  compared  to  the  baseline  implementation. 


EFFECT  OF  TARGET  NEUTRAL  ATOM  VELOCITY 

The  pre-collision  target  atom  velocity  distribution  has  a  great  impact  on  a  post-collision  CEX  ion  angular  distribution 
especially  for  elastic  collisions  that  involve  small  degree  of  momentum  transfer.  As  a  result,  the  differential  cross- 
section  at  larger  LAB  scattering  angles  is  smeared  out  and  approaches  the  pre-collision  target  atom  distribution. 
Figure  8  (a)  compares  the  LAB  differential  cross-section  obtained  for  different  pre-collision  target  atom  temperature. 
It  is  clearly  seen  that  the  CEX  population  with  large  angle  deflection  is  reduced  significantly.  However,  the  change 
in  the  differential  cross-section  is  only  attributed  to  slow  ions  having  post-collision  velocities  similar  to  their  initial 
velocities.  Since  the  slow  ions  are  easily  affected  by  the  electric  field,  the  change  in  the  differential  cross-section  does 
not  provide  any  information  to  what  degree  the  target  atom  velocity  distribution  affects  the  plume  simulation  results. 

In  order  to  confirm  that  the  stationary  target  gas  assumption  is  reasonable  in  plume  simulations,  the  speed  distri¬ 
bution  is  sampled  from  ion  population  after  a  single  elastic  collision  event.  Figure  8  (b)  shows  the  speed  distribution 
for  target  gas  temperatures  of  300  K  and  0  K.  For  almost  the  entire  range  of  the  kinetic  energy,  the  distribution  function 
is  nearly  the  same,  except  for  the  fluctuation  due  to  the  statistical  noise.  The  only  energy  range  that  is  significantly 
affected  by  the  target  gas  temperature  is  below  0.2  eV,  indicating  CEX  ion  trajectories  are  not  affected  by  the  assump¬ 
tion  as  long  as  there  is  a  potential  difference  between  a  cell  that  the  slow  ion  resides  in  and  neighboring  computational 
cells  that  is  more  than  0.2  V.  This  is  often  the  case  in  electric  propulsion  plumes,  especially  near  the  thruster  exit 
where  the  CEX  ions  are  mostly  generated. 


FIGURE  8.  (a)  LAB  differential  cross-sections.  The  theoretical  differential  cross-section  is  computed  for  station¬ 
ary  target  gas  (0  K)  whereas  the  one  sampled  from  particles  are  obtained  with  initial  target  particle  distribution  of 
Maxwellian  at  300  K.  (b)  Sampled  post-collision  speed  distribution  for  target  gas  temperature  of  300  K  and  0  K. 


CONCLUSION 


This  study  has  investigated  methods  to  accelerate  an  ion-atom  elastic  collision  calculation  that  includes  both 
momentum-  and  charge-exchange  processes  for  electric  propulsion  plume  simulations.  The  scattering  angles  are  pre¬ 
computed  through  a  classical  approach  with  ab  initio  spin-orbit  free  potential  and  are  stored  in  a  two-dimensional 
array  as  functions  of  impact  parameter  and  energy.  When  performing  a  collision  calculation  for  an  ion-atom  pair, 
the  scattering  angle  is  computed  by  a  table  lookup  and  multiple  linear  interpolations,  given  the  relative  energy  and 
randomly  determined  impact  parameter.  In  order  to  further  accelerate  the  calculations,  the  number  of  collision  cal¬ 
culations  is  reduced  by  properly  defining  an  effective  elastic  collision  cross-section,  <xe  1,  that  is  used  to  determine  if 
the  degree  of  momentum  exchange  is  insignificant.  The  effective  elastic  cross-section  is  a  function  of  energy  and  is 
used  to  obtain  the  maximum  impact  parameter  in  formulating  the  array  for  a  lookup  table.  In  the  MCC  method,  the 
target  atom  needs  to  be  sampled;  however,  it  is  confirmed  that  initial  target  atom  velocity  does  not  play  any  signifi¬ 
cant  role  in  typical  electric  propulsion  plume  simulations  such  that  the  sampling  process  is  unnecessary.  With  these 
implementations,  the  computational  run-time  to  perform  a  collision  calculation  is  reduced  by  71.4  percent  compared 
to  the  method  implemented  in  COLISEUM,  while  improving  accuracy  for  energies  that  the  previous  method  was  not 
optimized  for. 
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Xe++  Xe  Spin-Orbit  Free  Potential 


1 .  Fitting  potential  pair  into  Morse 
potential  form 

Vavg(r)  =  De(e2Hr‘~r)  -  2e^~r)) 


2.  Apply  statistical  weights 

V(r)  =  ^avg(r)Xvn^avg(r) 


[1]  F.  J.  Smith,  Physica  30,  497-504  (1964). 
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Normalized  CM  Deflection  Angle  (  x/n) 


Xe++Xe  Deflection  Function 
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Impact  Parameter,  A 


Deflection  angle  is  mostly 
affected  by  repulsive  potential 

•  Monotonically  decreasing  function 

•  Large  angle  at  short  range 

•  Small  angle  at  long  range 


CEX  Probability 

1  2  2  2 
Pcex(P)  =  ~ sin  A/7E+-sin  A rjn 


•  Highly  oscillatory  at  small  impact 
parameters  (PCEX=0.5) 

•  No  deflection  when  PCEX  has  some 
feature 
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Total  and  Differential  Cross-Section 


Introduction 


Numerical  Method 


Performance  Study 
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■k '  . 


•  Charge  exchange  cross-section1 

<jcex  =  87.3 -13.6  log(is) 

•  Effective  elastic  cross-section 

-  Indication  of  when  a  significant  momentum 
exchange  occurs 

-  Corresponds  to  a  cut-off  impact  parameter 


<jel  =  62.3 —  13.4  log(i?) 

•  Differential  cross-section 


91,  984-991  (2002). 


0  2  4  6  e 

Impact  Parameter,  A 


LAB  Scattering  Angle,  degree 
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Formulation  of  a  Lookup  Table 


Introduction 


Numerical  Method 


jfi 


Performance  Study 
Conclusion 
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•  Determine  scattering  angle  z(Er,b)  for  every  collision  pair 

-  Solving  classical  scattering  equation  is  slow 

-  Use  a  lookup  table  and  interpolate  to  collision  pair’s  relative 
energy  and  impact  parameter 


1 .  User  defined  parameters  Emin ,  Emax ,  NE ,  and  Nb 

2.  Get  max  impact  parameter  from  bma>(Er)  =  ^jEju 

3.  Solve  x  at  discrete  b  and  E  values 


E  —  E 

E  .=E  .  +  (  )  <  /  <  w- 1 

J  min  J  ^  5  J  E 

b,=i  bmax  ,  0  <  /  <  V  - 1 

Nb~  1 
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Using  a  Lookup  Table 
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Implementation  in  MCC  Method 
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Total  elastic  collision  cross-section  P‘ 


CEX 


<N.al  =  2(7CEX 
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Further  Improvement  is  Possible 
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Total  elastic  collision  cross-section 
Total  —  Tt.x  0-5cre| 

Charge  exchange  probability 
Pcex, i  —  0.5  and  PCEX j2  —  1  -0 
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Performance  Study 


•  Apply  collisions  to  1,146,700  particles  injected  at  270 
eV  (19,920  m/s) 

•  Repeat  calculations  for  100  times 

•  Only  reported  collision  calculation  time 

•  Baseline:  reduced  model  based  on  differential  cross- 
section  for  270  eV  LAB  energy  ions 


Model 

Wall  Time  (s) 

%  Baseline 

Baseline  (Sampling  from  curve-fit  to  differential  a) 

61.67 

100.0 

(1)  Lookup  Table 

45.75 

74.2 

(2)  Effective  Elastic  Collision  o  +  (1) 

36.21 

58.7 

(3)  Stationary  Background  Gas  4-  (2) 

17.63 

38.6 
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Effect  of  Target  Neutral  Atom  Velocity 
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Does  the  assumption  made  to  pre-collision  neutral 
velocity  have  any  significant  effect  to  the  plume 
simulation? 


•  Differential  cross-section  is  altered  significantly 

-  Especially  the  large  angle  scattered  population 

•  This  population  only  has  small  energy  (<0.2eV) 
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Introduction 
Numerical  Method 
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Conclusion 


Conclusion 


•  Acceleration  of  elastic  collision 
calculation  without  losing  accuracy 
from  high-fidelity  model 

-  Effective  elastic  collision  cross-section 

-  Formulation  of  a  lookup  table 

•  Uniformly  distributed  E  and  b 

-  Implementation  with  based  on  three 

different  cross-sections 

•  Stationary  background  gas 
assumption  is  sufficient  for  plume 
simulations 


•  >70%  improvement  in  speed 
compared  to  a  reduced  model 
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Thank  you 
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